Houston Storm Blackouts

Identifying the impacts of extreme weather

Skill Demonstration: Manipulating vector & raster data | Spatial joins across data types | Spatial analysis across data types

Photo source: soa.org

Purpose

This repository demonstrates my analysis of the 2021 blackouts in Houston resulting from an extreme winter storm.

Skills and Techniques Demonstrated

• Manipulating raster data
• Manipulating vector data
• Spatial joins across data types
• Spatial analysis across data types

# libraries
library(terra)
library(sf)
library(tidyverse)
library(here)
library(tmap)
library(tmaptools)
library(viridisLite)
library(stars)

# starting with a bounding box of the houston area to minimize processing
bbox <- st_bbox(c(xmin = -96.5,
                  ymin = 29,
                  xmax = -94.5,
                  ymax = 30.5))

# load data
## open street map buildings
osm_roads <- st_read(here("data", "gis_osm_roads_free_1.gpkg"), 
                 query = "SELECT * FROM gis_osm_roads_free_1 WHERE fclass='motorway'")

osm_buildings <- st_read(here("data", "gis_osm_buildings_a_free_1.gpkg"), 
query = "SELECT * FROM gis_osm_buildings_a_free_1 WHERE (type IS NULL AND name IS NULL) OR type in ('residential', 'apartments', 'house', 'static_caravan', 'detached')")

# check / convert buildings CRS

if(crs(osm_buildings) == "EPSG: 3083"){
  print("CRS correct")
} else{
  warning("Updated CRS to EPSG 3083")
  osm_buildings <- st_transform(osm_buildings, "EPSG: 3083")
}

# check / convert roads CRS
if(crs(osm_roads) == "EPSG: 3083"){
  print("CRS correct")
} else{
  warning("Updated CRS to EPSG 3083")
  osm_roads <- st_transform(osm_roads, "EPSG: 3083")
}

# light data
light_02_07a <- rast(here("data", "VNP46A1", "VNP46A1.A2021038.h08v05.001.2021039064328.tif"))

light_02_07b <- rast(here("data", "VNP46A1", "VNP46A1.A2021038.h08v06.001.2021039064329.tif"))

light_02_16a <- rast(here("data", "VNP46A1", "VNP46A1.A2021047.h08v05.001.2021048091106.tif"))

light_02_16b <- rast(here("data", "VNP46A1", "VNP46A1.A2021047.h08v06.001.2021048091105.tif"))


# load census tract geometries
census_tracts <- st_read(here("data", "ACS_2019_5YR_TRACT_48_TEXAS.gdb"),
                  layer = "ACS_2019_5YR_TRACT_48_TEXAS")

# check / convert census tract data CRS
if(crs(census_tracts) == "EPSG: 3083"){
  print("CRS correct")
} else{
  warning("Updated CRS to EPSG 3083")
  census_tracts <- st_transform(census_tracts, "EPSG: 3083")
}

# load income data
income <- st_read(here("data", "ACS_2019_5YR_TRACT_48_TEXAS.gdb"),
                  layer = "X19_INCOME")

# join the data tables
census_income <- left_join( census_tracts, income, 
                            join_by(GEOID_Data == GEOID))

# check / convert census_income data CRS
if(crs(census_income) == "EPSG: 3083"){
  print("CRS correct")
} else{
  warning("Updated CRS to EPSG 3083")
  census_income <- st_transform(census_income, "EPSG: 3083")
}

Part 1: Mapping the Houston Blackouts

# merge the 02-07 light rasters
light_02_07 <- merge(light_02_07a, light_02_07b)

# merge the 02-16 light rasters
light_02_16 <- merge(light_02_16a, light_02_16b)

Maps of Houston nightlight intensity before and after the storm.

#crop the night light data to the Houston area
h_light_0207 <- terra::crop(light_02_07, bbox)
summary(h_light_0207)

h_light_0216 <- terra::crop(light_02_16, bbox)

tmap_mode("plot")

# February 7, 2021
map_0207 <- tm_shape(h_light_0207)+
  tm_raster(breaks = c(0, 0.2, 1, 3, 5, 10, 100, 200, 10000, 100000),
    palette = magma(10),
            title = "Night Light Intensity")+
  tm_layout(legend.position = c("left", "top"),
            legend.text.size = 0.6,
            legend.text.color = "white", 
            legend.title.size = 0.8, 
            legend.title.color = "white",
            legend.bg.color = "black",       #
            legend.bg.alpha = 0.7,  
            outer.bg.color = "transparent",
            main.title = "February 07, 2021")+
  tm_compass(position = c("right", "bottom"))+
  tm_scale_bar(position = c("left", "bottom"))

tmap_save(map_0207, here("outputs", "map_0207.png"), bg = "transparent", dpi = 400)

# Feb 16, 2021
map_0216 <- tm_shape(h_light_0216)+
  tm_raster(breaks = c(0, 0.2, 1, 3, 5, 10, 100, 200, 10000, 100000),
    palette = magma(10),
            title = "Night Light Intensity")+
  tm_layout(legend.position = c("left", "top"),
            legend.text.size = 0.6,
            legend.text.color = "white", 
            legend.title.size = 0.8, 
            legend.title.color = "white",
            legend.bg.color = "black",       #
            legend.bg.alpha = 0.7,  
            outer.bg.color = "transparent",
            bground.color = "transparent",
            main.title = "February 16, 2021")+
  tm_compass(position = c("right", "bottom"))+
  tm_scale_bar(position = c("left", "bottom"))

tmap_save(map_0216, here("outputs", "map_0216.png"), bg = "transparent", dpi = 400)

Maps Interpretation

Note the change in the urban center. The general pattern of blackouts are visible in the form of orange encroaching into the yellow area. This is apparent despite the overall increase in light intensity across the area.

Part 2: Identifying impacted homes and census tracts based on changes in night light intensity.

Homes impacted by blackout

# use h_light_0207 to make a Houston area vector for later use
houston <- st_as_stars(h_light_0207)

# vectorize
houston <- st_as_sf(houston)

# check / convert Houston cookie cutter data CRS

if(crs(houston) == "EPSG: 3083"){
  print("CRS correct")
} else{
  warning("Updated CRS to EPSG 3083")
  houston <- st_transform(houston, "EPSG: 3083")
}

# dissolve the grid into cookie cutter
houston <- st_union(houston)

# Calculate light intensity difference
light_diff <- light_02_07 - light_02_16
summary(light_diff)

# Calculate light intensity difference

## define reclass matrix
rcl <- matrix(c(-20900, -200, 1,# group 1 ranges from -20900 (min) to -200
                -201, 3775, NA), # group 2 ranges from -200 to 3776 (max)
                ncol = 3, byrow = TRUE)

# use reclass matrix to reclassify light_diff raster
blackouts <- classify(light_diff, rcl = rcl)

# change reclasssed values into factors
values(blackouts) <- as.factor(values(blackouts))
summary(blackouts)

# Crop blackout raster to Houston area
blackouts_cropped <- terra::crop(blackouts, bbox)
summary(blackouts_cropped)


# Vectorize light raster
## first, convert it to a stars object
blackouts_stars <- st_as_stars(blackouts_cropped)

# have a look at it
summary(blackouts_stars)

# vectorize!!!
blackouts_vector <- st_as_sf(blackouts_stars)

# check / convert blackout vector CRS
if(crs(blackouts_vector) == "EPSG: 3083"){
  print("CRS correct")
} else{
  warning("Updated CRS to EPSG 3083")
  blackouts_vector <-  st_transform(blackouts_vector, "EPSG: 3083")
}

# load county data as a basemap
counties <- st_read(here("data", "tl_2023_us_county", "tl_2023_us_county.shp"))

# check / convert county vector CRS
if(crs(counties) == "EPSG: 3083"){
  print("CRS correct")
} else{
  warning("Updated CRS to EPSG 3083")
  counties <-  st_transform(counties, "EPSG: 3083")
}

# crop the county data to the Houston area
houston_counties <- st_filter(x = counties, y = houston,
                               .predicate = st_intersects)

# Make a 200 m buffer around roads
## check for meters as units
st_crs(osm_roads)$units 

# buffer the roads by 200m
highway_buffer <- st_buffer(osm_roads, dist = 200)

# combine the buffer polygons
highway_buffer_union <- st_union(highway_buffer)

## Exclude highway buffer from blackout mask
blackouts_mask_no_highways <- st_difference(blackouts_vector, highway_buffer_union)

# Spatial geometry filtering to count homes impacted by blackouts
blackout_homes <-  st_filter(x = osm_buildings, y = blackouts_mask_no_highways, .predicate = st_intersects)

bbox <-  st_bbox(blackout_homes)

map_blackout_homes <- tm_shape(counties, bbox = bbox)+
  tm_polygons(col = "#1e0057",
              fill = "#1e0057",
              fill_alpha = 0.7)+
  tm_shape(houston_counties) +
  tm_polygons(col = "#1e0057",
              fill = "#3b00ab") +
  tm_shape(highway_buffer_union)+
  tm_polygons(col = "#1e0057",
              fill = "#1e0057",
              fill_alpha = 0.5) +
  tm_shape(blackout_homes)+
  tm_polygons(col = "#fd007b")+
  tm_layout(outer.bg.color = "transparent",
            bg.color = "#000327")

tmap_save(map_blackout_homes, here("outputs", "map_blackout_homes.png"), dpi = 400)

Census tracts impacted by blackout

## Crop census data to the Houston area
census_income <- st_filter(x = census_income, y = houston,
                           .predicate = st_intersects)

# filter down to the affected tracts
blackout_tracts <- st_filter(x = census_income, y = blackouts_mask_no_highways,
                             .predicate = st_intersects)

# filter down to the unaffected tracts 
## (st_difference isn't working so this is my work-around)

# find the indices of blackout polygons
blackout_indices <- st_intersects(census_income, blackouts_mask_no_highways)

# filter by those indices
unaffected_tracts <- census_income[!sapply(blackout_indices, length) > 0, ]

bbox2 <- st_bbox(blackout_tracts)

map_blackouts_tracts <- tm_shape(counties, bbox = bbox2)+
  tm_polygons(col = "#1e0057",
              fill = "#1e0057",
              fill_alpha = 0.7)+
  tm_shape(houston_counties) +
  tm_polygons(col = "#1e0057",
              fill = "#3b00ab") +
  tm_shape(highway_buffer_union)+
  tm_polygons(col = "#1e0057",
              fill = "#1e0057",
              fill_alpha = 0.5) +
  tm_shape(blackout_tracts)+
  tm_polygons(col = "#fd007b",
              fill = "#fd007b",
              fill_alpha = 0.5)+
  tm_layout(outer.bg.color = "transparent",
            bg.color = "#000327")

tmap_save(map_blackouts_tracts, here("outputs", "map_blackouts_tracts.png"), dpi = 400)

Homes likely to have experiencced the blackout are depicted as pink polygons in the above figure. An estimated 81914 homes experienced blackouts.

Here the census tracts likely to have experienced the blackout are visible as pink polygons overlaid on the Houston area.

Part 3: Comparing the median household income among impacted and un-impacted census tracts.

# Add a grouping variable to your data beforehand if needed
blackout_tracts$status <- "Impacted"
unaffected_tracts$status <- "Unaffected"

# Combine the data
combined_tracts <- rbind(blackout_tracts, unaffected_tracts)

plot_blackout_tracts <- ggplot(data = combined_tracts, aes(x = B19013e1, fill = status)) +
  geom_histogram(data = blackout_tracts, alpha = 0.5, bins = 30) +
  geom_histogram(data = unaffected_tracts, alpha = 0.5, bins = 30) +
  scale_fill_manual(
    name = "Census Tract Status",  # Legend title
    values = c("Impacted" = "#fd007b", "Unaffected" = "#3b00ab"),
    labels = c("Impacted by Blackout", "Not Impacted")) +  # Legend labels
  labs(x = "Median Household Income ($)",
       y = "",
       title = "Distribution of Median Incomes for Census Tracts Impacted by Blackouts") +
  theme_bw() +
  theme(panel.background = element_rect(fill = 'transparent'),
        plot.background = element_rect(fill = 'transparent', color = NA))

ggsave(here("outputs", "plot_blackout_tracts.png"), plot_blackout_tracts,
       height = 45, width = 7, dpi = 400, bg = "transparent")

Figure Interpretation

Above are depicted wide spread of incomes impacted by the blackouts. Note that lower median household incomes are more highly represented than un-impacted tracts.

Discussion

In the above analysis and accompanying figures, the impacts of the storm are made apparent as a visible change in night light intensity in the hoston metro area. An estimated 81914 homes experienced blackouts according to this analysis. Further, these results suggest that blackouts occured independently of income level, as indicated by the median household incomes across census blocks. This is somewhat surprising as one might expect households with higher income to overcome barriers and power their homes (and lights). This trend implies the pervasive nature of the blackouts without readily available alternate power sources.

Acknowledgements

This analysis and workflow was originally created by Dr. Ruth Oliver of the Bren School of Environmental Science & Management for the Masters of Environmental Data Science course, Geospatial Analysis and Remote Sensing. My initial work through of this analysis was conducted as part of a homework assignment for this class, and I later polished up this repository.

Data

The data used in this analysis is too large to be hosted on GitHub. Instead, download the data here as a zipped folder, unzip, and move into the R project manually.

Citations

Data	Link
Night Lights	https://ladsweb.modaps.eosdis.nasa.gov/
Open Street Map Roads	https://planet.openstreetmap.org/
Open Street Map Buildings	https://planet.openstreetmap.org/
Socioeconomic Data	https://www.census.gov/programs-surveys/acs